""" 点密度地图表达了给定区域内点状符号的密度 """
import shapefile
import random
import pngcanvas


# 先判断是否在多边形内部
def point_in_poly(x, y, poly):
    # 判断是否是顶点
    if (x, y) in poly:
        return True
    # 判断是否在多边形边上
    for i in range(len(poly)):
        p1 = None
        p2 = None
        if i == 0:
            p1 = poly[0]
            p2 = poly[1]
        else:
            p1 = poly[i-1]
            p2 = poly[i]
        if p1[1] == p2[1] and p1[1] == y and min(p1[0], p2[0]) < x < max(p1[0], p2[0]):
            return True
    n = len(poly)
    inside = False
    p1x, p1y = poly[0]
    for i in range(n+1):
        p2x, p2y = poly[i % n]
        if y > min(p1y, p2y):
            if y <= max(p1y, p2y):
                if x <= max(p1x, p2x):
                    if p1y != p2y:
                        xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
                        if p1x == p2x or x <= xints:
                            inside = not inside
            p1x, p1y = p2x, p2y
    if inside:
        return True
    else:
        return False


# 地理坐标转换为屏幕坐标
def world2screen(bbox, w, h, x, y):
    minx, miny, maxx, maxy = bbox
    xdist = maxx - minx
    ydist = maxy - miny
    xratio = w/xdist
    yratio = h/ydist
    px = int(w - ((maxx - x) * xratio))
    py = int((maxy - y) * yratio)  # 为了适应pngcanvas的坐标系，对y做了翻转
    return px, py


# 主程序
def main():
    # 打开人口普查shp文件
    r = shapefile.Reader("shapefiles/GIS_CensusTract/GIS_CensusTract_poly")
    # print(r)
    # 45 shapes(type 'POLYGON')
    # 45 records(24 fields)

    # 获取人口记录索引
    pop_index = None
    dots = []
    for i, f in enumerate(r.fields):
        if f[0] == 'POPULAT11':
            # 声明删除标记
            pop_index = i - 1
    # 计算点密度并绘制相关点
    for sr in r.shapeRecords():
        population = sr.record[pop_index]
        # 密度比例 1 个点代表 100 人
        density = population/100
        found = 0
        # 随机绘制点，直到密度达到制定比率
        while found < density:
            minx, miny, maxx, maxy = sr.shape.bbox
            x = random.uniform(minx, maxx)
            y = random.uniform(miny, maxy)
            if point_in_poly(x, y, sr.shape.points):
                dots.append((x, y))
                found += 1

    # 为输出PNG图片做准备
    # 设置输出图片尺寸
    width = 600
    height = 400
    c = pngcanvas.PNGCanvas(width, height)
    c.color = (255, 0, 0, 0xff)
    for d in dots:
        # 将使用*d标记扩展(x,y)元祖
        x, y = world2screen(r.bbox, width, height, *d)
        c.filled_rectangle(x-1, y-1, x+1, y+1)
    # 绘制人口普查区域
    c.color = (0, 0, 0, 0xff)
    for s in r.iterShapes():
        pixels = []
        for p in s.points:   # for x,y in s.points
            pixel = world2screen(r.bbox, width, height, *p)
            pixels.append(pixel)
        c.polyline(pixels)

    # 保存图片
    img = open("DotDensity.png", "wb")
    img.write(c.dump())
    img.close()
    r.close()


if __name__ == '__main__':
    main()
